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ABSTRACT 

Fluctuations in a stellar system’s gravitational field cause the orbits of stars to evolve. The resulting 
evolution of the system can be computed with the orbit-averaged Fokker-Planck equation once the 
diffusion tensor is known. We present the formalism that enables one to compute the diffusion tensor 
from a given source of noise in the gravitational field when the system’s dynamical response to that 
noise is included. In the case of a cool stellar disc we are able to reduce the computation of the 
diffusion tensor to a one-dimensional integral. We implement this formula for a tapered Mestel disc 
that is exposed to shot noise and find that we are able to explain analytically the principal features of 
a numerical simulation of such a disc. In particular the formation of narrow ridges of enhanced density 
in action space is recovered. As the disc’s value of Toomre’s Q is reduced and the disc becomes more 
responsive, there is a transition from a regime of heating in the inner regions of the disc through the 
inner Lindblad resonance to one of radial migration of near-circular orbits via the corotation resonance 
in the intermediate regions of the disc. The formalism developed here provides the ideal framework 
in which to study the long-term evolution of all kinds of stellar discs. 

Subject headings: Galaxies, dynamics, evolution, diffusion 


1. INTRODUCTION 

Many, perhaps all, stars are born in a stellar disc. Ma¬ 
jor mergers destroyed some discs quite early in the his¬ 
tory of the Universe, but many others have survived to 
the present day, including the disc of which the Sun is 
a part. Hence an understanding of the dynamics and 
evolution of stellar discs is an essential ingredient of cos¬ 
mology. Conversely, cosmology provides the framework 
within which disc dynamics should be studied because 
dark-matter halos make large contributions to the gravi¬ 
tational fields in which discs move, and dark-matter sub¬ 
structures are major contributors to the gravitational 
noise to which discs are exposed. 

Serious study of disc dynamics got underway in the 
1960s with seminal works by Lin, Shu, Goldreich, Toornre 
and Lynden-Bell. Although some important insights 
were gained at that stage, fundamental questions were 
left open. While the earliest work was almost entirely 
analytic in nature, numerical simulations of stellar discs 
became more important over time, and revealed impor¬ 
tant aspects of disc dynamics that were hard to under¬ 
stand analytically. In particular it emerged that discs 
that are completely stable at a linear level nevertheless 
develop spiral structure that eventually grows to ampli¬ 
tudes of order unity, so the disc becomes something like 
a barred spiral galaxy (Sellwood 2012, hereafter S12). 
Sellwood & Carlberg (2014) recently offered a convinc¬ 
ing explanation of this phenomenon that hinges on the 
fact that resonances localise the impact that a fluctua¬ 
tion has on a disc. This localisation is the major focus 
of this paper. 

Self-gravitating discs are responsive dynamical sys- 
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terns, in which (a) rotation provides an abundant sup¬ 
ply of free energy, and (b) resonances play a key role. 
The ready availability of free energy leads to some stim¬ 
uli being powerfully amplified, while resonances localise 
the dissipation of free energy with the result that even 
a very small stimulus can result in a disc evolving to an 
equilibrium that is materially different from the one from 
which it started. 

The stimuli to which discs respond are various sources 
of gravitational noise. They include, Poisson noise aris¬ 
ing from the finite number of stars in a disc, Poisson 
noise arising from the finite number of giant molecular 
clouds in the interstellar medium, and Poisson noise aris¬ 
ing from the finite number of massive sub-halos around 
a galaxy. Spiral arms in the distribution of gas provide 
another source of gravitational noise, while the rotating 
gravitational field of a central bar constitutes a source of 
stimulus that is more systematic than noisy. The history 
of a real stellar disc will largely comprise responses to all 
these stimuli. 

In the solar neighbourhood at least three distinct man¬ 
ifestations of such responses are evident: 

(i) The random velocity of each coeval cohort of 
stars increases with the cohort’s age (Wielen 1977; 
Aumer & Binney 2009). 

(ii) The velocity distribution at the Sun contains sev¬ 
eral “streams” of stars (Delmen 1998). Each 
such stream contains stars of various ages and 
chemistries that are all responding to some stim¬ 
ulus in a similar way (Famaey et al. 2005). 

(iii) In the two-dimensional space in which one coordi¬ 
nate is angular momentum and the other is a 
measure of a star’s radial excursions, such as the 
radial action J r , the density of stars shows elon¬ 
gated features. The density of stars is depressed 
near J r = 0 but enhanced at larger J r in such 
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a way that the whole region of disturbed stellar 
density forms a curve that is consistent with be¬ 
ing a curve on which a resonant condition such as 
2f — = constant holds (Sellwood 2010; McMil¬ 

lan 2013). We shall call a feature of this type a res¬ 
onance ridge. Sellwood & Carlberg (2014) have 
argued that resonance ridges play a crucial role in 
the long-term dynamics of stellar discs. 

Numerical simulations of stellar discs are extremely 
challenging because the near two-dimensional geometry 
of discs combined with their responsive nature causes dis¬ 
creteness noise to be dynamically important unless the 
number of particles employed exceeds ~ 200 000. Hence 
only recently has it become straightforward to simulate a 
disc with a sufficient number of particles for discreteness 
noise to be dynamically unimportant for many dynamical 
times (S12). It is particularly hard to simulate accurately 
a disc that is embedded in a cosmological simulation and 
thus exposed to cosmic noise. Moreover, the utility of 
a simulation is greatly increased if one understands an¬ 
alytically why it evolves the way it does. A goal of this 
paper is to show the extent to which perturbation theory 
explains a phenomenon - resonance ridges - that is seen 
in both numerical simulations and surveys of the solar 
neighbourhood. 

Perturbation theory is much more than a device for 
computing approximate solutions to equations: through¬ 
out physics it provides the conceptual framework we 
use to understand phenomena. Examples include the 
concepts of a free particle and an interaction in parti¬ 
cle physics, a phonon and a gravity wave in condensed- 
matter physics, semi-major axis and eccentricity in plan¬ 
etary dynamics, and so on. The natural way to increase 
our understanding of the dynamics of stellar discs is to 
practise the application of perturbation theory to these 
systems, so we may gain insight into how these fascinat¬ 
ing systems work, and learn how one can think about 
them most profitably. 

Kalnajs (1971) laid the foundations of perturbation 
theory for stellar discs. The theory is based on the use 
of angle-action coordinates - the coordinates that were 
introduced to understand the dynamics of the solar sys¬ 
tem. These coordinates are being increasingly used to 
build equilibrium models of hot and cold stellar systems 
(Binney 2010, 2014; Piffl et al. 2014), and to study the 
dynamics of star streams (Helmi & White 1999; Sellwood 
2010; McMillan 2013; Eyre & Binney 2011; Sanders & 
Binney 2013). Binney & Lacey (1988) used these coor¬ 
dinates to derive the orbit-averaged Fokker-Planck equa¬ 
tion for a stellar disc. However, they did not consider 
the origin of the fluctuations in the gravitational poten¬ 
tial that drive stellar diffusion. Weinberg (2001a) di¬ 
vided the driving fluctuations into the contribution from 
some external stimulus, and the self-consistent dynami¬ 
cal response of the system itself to the stimulus. Wein¬ 
berg’s treatment was adapted to systems that are spher¬ 
ical when unperturbed, while here we restrict ourselves 
to razor-thin discs, in which case the construction of the 
angle-action coordinates is trivial. 

The paper is organised as follows. Section 2 recalls 
from Binney & Lacey (1988) and Weinberg (2001a) 
the general principles of secular evolution, the orbit- 
averaged Fokker-Planck equation, and the use of a set 


of biortlionornral potential-density pairs to compute the 
diffusion tensor that is jointly generated by an exter¬ 
nal stimulus and the system’s response to this stimulus. 
Section 3 specialises this formalism to a razor-thin, cool 
disc by introducing a set of basis functions that comprise 
localised, tightly-wound spirals. Using these basis func¬ 
tions we are able to reduce the computation of the diffu¬ 
sion tensor, which in principle requires a double sum over 
basis functions, to a single integral over radial wavenum¬ 
bers. In Section 4 we compute the diffusion tensor for 
a tapered Mestel disc that is excited by shot noise, and 
show that the resulting predictions for the disc’s evolu¬ 
tion reproduce the main features of the N-body simula¬ 
tions reported by S12. Finally, we conclude in Section 5. 

2. FLUCTUATIONS AND SECULAR EVOLUTION 

To zeroth order, stellar discs are systems that have 
achieved statistical equilibrium within an axisymmetric 
gravitational field that arises not only from their mass 
but also from mass contained in other components of 
the galaxy, especially the bulge and the dark halo. The 
Hamiltonian associated with the field is to a good approx¬ 
imation integrable, so all orbits may be assumed to admit 
three isolating integrals, which we take to be the actions: 
J<t> = L z is the angular momentum about the field’s sym¬ 
metry axis; J r , which quantifies the amplitude of a star’s 
radial oscillations; and J 2 , which quantifies oscillations 
perpendicular to the field’s equatorial plane (Born 1960; 
Binney & Tremaine 2008). On account of the integra- 
bility of the gravitational field and Jeans’ theorem, we 
can assume that at each instant the disc’s distribution 
function (DF) is a function /(J,i) of the actions only, 
rather than a general function on phase space /(J, 9 , f), 
which has dependence on the variables that are canoni¬ 
cally conjugate to the actions, namely the angle variables 

Any fluctuation in the gravitational field causes each 
star to deviate from its original orbit J and to settle 
after the fluctuation has died away on another orbit 
J' = J + A. Hence fluctuations cause stars to diffuse 
through action space. Since initially action space is pop¬ 
ulated by stars only along the axis, this diffusion raises 
the density of stars away from this axis by populating 
orbits with distinctly non-zero J r . As a consequence, 
the velocity dispersion within the disc rises, so fluctu¬ 
ations “heat” the disc. Stars also diffuse along the 
axis. Since such diffusion merely transfers stars from one 
nearly circular orbit to another of a different radius, this 
component of diffusion is called radial migration. Radial 
migration does not heat the disc, and given that the den¬ 
sity of stars does not vary rapidly along the J$ axis, it can 
easily go unnoticed. However, chemical evolution within 
the disc establishes a radial gradient in metallicity, and 
radial migration is most readily detected through its in¬ 
teraction with this gradient (Sellwood & Binney 2002): 
radial migration tends to erase the correlation between 
the ages and metallicities of stars near the Sun by bring¬ 
ing to the Sun both old, nretal-rich stars formed at small 
radii and young metal-poor stars formed at large radii. 

Fundamentally fluctuations drive the long term (“sec¬ 
ular”) evolution of discs in much the same way that 
they drive the much better understood secular evolution 
of globular clusters, but resonances are unimportant in 
globular clusters and dominant in discs. As indicated in 
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the Introduction, resonances localise the impact of fluc¬ 
tuations and give rise to ridges in action space that are 
the primary focus of this paper. However, Sellwood & 
Binney (2002) pointed out that at the corotation reso¬ 
nance, stars are scattered at constant radial action, i.e., 
parallel to the axis. So scattering at corotation does 
not give rise to a ridge and is likely to go overlooked if 
one does not pay attention to the nretallicities of stars. 
While at corotation only J$ changes, at a Lindblad reso¬ 
nance both J r and J$ change, so radial migration is not 
confined to corotation, as is often stated. 


By the Wiener-Khinchin theorem, c m (J, v) is the power 
spectrum of the stationary random variable 
Equation (7) tells us that diffusion is driven by reso¬ 
nances because it implies that the rate at which stars 
diffuse from action J is proportional to the power that 
the fluctuating field has at any of the orbit’s character¬ 
istic frequencies m • fi(J). Hence, if the fluctuations are 
confined to a narrow frequency range, perhaps because 
they are associated with spiral arms, stars that respond 
strongly to them will be located at only a few points in 
action space. 


2.1. Orbit-averaged Fokker-Planck equation 


Since we are imagining that stars are conserved, the 
equation that governs the secular evolution of the DF 
takes the form 


dt dJ ’ 


(1) 


where F is the diffusive flux of stars in action space. If 
P(J, A) is the rate of increase with time of the probabil¬ 
ity that a star scatters from J to J + A, then F is given 
by (Binney & Lacey 1988) 


Ft = /A,, 



( 2 ) 


where the first- and second-order diffusion coefficients are 


Aj(J) = /d 3 A AiP(J, A) 
A£(J) = J d 3 A A,AjP(J, A). 


(3) 


Binney & Lacey (1988) showed that in the relevant cir¬ 
cumstances the first- and second-order diffusion coeffi¬ 
cients are related by 


A, = 


1 d _^l 

2 dJj ’ 


( 4 ) 


so the diffusive flux can be written entirely in terms of 
the second-order coefficients: 


2.2. A basis-function expansion 

We will find it expedient to expand the fluctuating po¬ 
tential ip(x,t) in a set of basis functions, the members of 
which are enumerated by an index q: 

^(x, t)=^2 b q (t)ip iq) (x), (9) 

4 

where b q (t) is a random variable. Following Kalnajs 
(1971), we require our basis potentials to be orthonor¬ 
mal to the densities that generate them, so we 

have 

V 2 V> (p) = 4 t xGpW and f d 3 xp^(x)[^ (9) (x)]* = -S pq . 

( 10 ) 

Now we have 

V’mGM) = ■jywj f d 3 0e -im 'V(0, J ,t) 

71 J ( 11 ) 

= XlMOV’mG 1 ). 

P 

where 

V^(J) = j d 3 0e- im 'V p) [x(0,J)]. (12) 

Hence the required power spectrum is 

= 5> P9 (^ } ( J)V& ] *(J), (13) 

pq 


Fi 


= -i^r.2L 
2 '’dJ* 


(5) 


By expanding ^(x, t), the fluctuating part of the grav¬ 
itational potential, in angle-action variables, 

^(x,t) = ip(0,3,t) = ^V’m(J,f)e lm ' 0 , (6) 

m 


where 

/ OO 

d t e WT b p (t)b*(t - t) (14) 

-OO 

is the Fourier transform of the cross-correlation of the 
amplitudes of the p and q basis functions. 

Below we shall require an expression for B pq (o) in 
terms of the Fourier transform 


Binney & Lacey (1988) showed that the second-order dif¬ 
fusion coefficients are related to the fluctuations in the 
potential by 

Km = y m»TOjC m (J, m • »(J)), (7) 

m 

where an overline indicates an ensemble average and c m 
is the Fourier transform with respect to time of the auto¬ 
correlation of the m Fourier component of the potential 
/>00 

Cm(J,i/)= / dr e 1UT (J, (J, t — r). (8) 


b p (v) = J d te wt b p (t) (15) 

of b p {t). If ip(t) is a stationary random process, then it 
is straightforward to show that 

bp(y)b* q {v') = 2t:5(v - iy')B pq (iy). (16) 

2.3. Bare and dressed stimuli 

As we indicated in the Introduction, a stellar disc is ex¬ 
posed to several sources of fluctuations. The issue that 
we now have to confront is that the fluctuation if in the 


— OO 
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potential that stars experience, which is what appears 
in the above formulae, differs from the original stimu¬ 
lation, i/j e , because the disc has non-negligible mass, so 
through Poisson’s equation it makes a contribution i/> s to 
the actual gravitational potential ip (Weinberg 2001a). 
We shall refer to ip e as the “bare” stimulus and to 

m = r(t)+r(t) (i7) 

as the “dressed” stimulus. We now seek a relationship 
between the dressed and bare stimuli. 

Let ip' be the change in the potential that the disc 
would generate if its particles moved in the sum of the 
unperturbed potential and the stimulating potential ip e . 
Then t//(x) is linearly related to V’ e ( x ) so f° r each time- 
lapse r there is a linear response operator M(r) that 
connects these functions 


equation (23) by f d te lvt yields 

b r {v) = a r {v) + ^2 M rq {v)b q {v). (25) 

Q 

Hence 

b(z/) = [I - M(i/)] _1 a(^), (26) 

where boldface implies vectors and matrices indexed with 
13- 

Equations (16) and (26) enable us to relate B. pq to the 
basis coefficients of the stimulating field 

B (y) = — I d v' [I—M(^)] _ 1 a(^) ® a*(i/')[I-M t (i/ , )] _1 - 
2 ttJ 

(27) 

We will show below that analogously to equation (16) 


ip'(t) = ( d t'M(t — t')ijj e (t'). (18) 

J — OO 

Since the mass of the disc actually contributes to the 
potential in which its particles move, changes in the disc’s 
potential at a early time t' contribute alongside ip e (t') to 
the disturbance of disc particles at later times, so the 
fluctuating component of the disc’s potential ip s satisfies 

ip s {t) = I dt'M(t - t')[ip s (t') + ip e (t’)}. (19) 

J — OO 

Inserting this expression into equation (17), we obtain 
ip(t) = ip e (t) + f d t'M(t — t')ip(t'). (20) 

J —OO 

In this equation the potentials are functions of x as well 
as t and — is an operator that maps one function 
of space onto another. The basis functions introduced 
above reduce this operator to a matrix, so when we write 


V> e CM) = ^2a P (t)ip ip \x.) 

P 

M) = 5Z 6 p(^ ( p) ( x )’ 


equation (20) can be written 

X]M^ (p) ( x ) [ a «W (9) ( x ) 

P Q 

+ [ dt'M(t — t')b q (t , )ip ( - q \x.) . 

J —OO J 

( 22 ) 

We multiply both sides of this equation by — / d 3 x [p^]* 
and with equation (10) obtain 

b r (t) = a r (t)+ ^2 [ d t'M rq (t — t')b q (t'), (23) 

q J-°° 

where 

M rq (t -t') = - J d 3 x [p^ (x)]* M(t - t’) V> (9) (x). (24) 


The temporal convolution can be reduced to a multipli¬ 
cation by taking a Fourier transform: multiplication of 


a P (v)a*(v') = 2tt5(v - v’)A pq (v). (28) 

Hence equation (13) can be written 

On (J, I /)=X!^ ) ( J )^ ) *( J ) 

pq 

x {[I - M(^)] _1 A(^)[I - M^I/)]- 1 } . 

I J pq 

(29) 

Our derivation of the dressed secular diffusion coefficients 
sketched previously is based on the master equation (Bin- 
ney & Tremaine 2008) which led to the first- and second- 
order diffusion coefficients from equation (3). One can 
also recover these diffusion coefficients via a timescale de¬ 
coupling of the collisionless Boltzmann equation (Wein¬ 
berg 2001a; Pichon & Aubert 2006; Chavanis 2012; Fou¬ 
vry et al. 2014a). Various sources of external pertur¬ 
bations can then be considered to induce secular evolu¬ 
tion (Weinberg 2001b; Aubert & Pichon 2007). 

3. APPLICATION TO A COOL, THIN DISC 

The simplest non-trivial context in which the above 
principles can be illustrated is the case of a cool razor- 
thin disc, i.e., a disc in which every star is confined to a 
plane by J z = 0 and orbits have only moderate eccentric¬ 
ities. In this case each unperturbed orbit is characterised 
by two numbers (J r , J^) or a point J in two-dimensional 
action space. The angular momentum is as ever a 
trivial function of (x, v), and since orbits have only mod¬ 
erate eccentricities, the epicycle approximation provides 
an adequate expression for J r (x, v). If k( Jq) denotes the 
radial epicycle frequency, the disc’s unperturbed DF can 
be taken to be an exponential exp (—nJ r /a^) in J r times 
a function of Jp that essentially controls the disc’s radial 
surface-density profile. Then within the epicycle approx¬ 
imation the velocity distribution at any point in the disc 
is a biaxial Gaussian with radial dispersion ay. Because 
gas tends to flow on nearly closed orbits, stars are born 
on nearly circular orbits, i.e., along the J$ axis of action 
space, and diffuse from there in to the body of action 
space. We shall show that on account of resonances, this 
diffusion can form ridges in action space. 

3.1. Choice of the basis 

Since we are working in two dimensions, the basis func¬ 
tions (x) become functions ip^ (R, <j>) of plane polar 
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coordinates that are orthonormal to the generating sur¬ 
face densities E ^(R, <j>). Our problem is simplified if we 
can choose the basis pA' 1 such that the response operator 
M is diagonal and we now show that this is possible. 

It is well known that the potential generated via Pois¬ 
son’s equation by a tightly wound spiral wave is itself a 
spiral wave with the same wavevector k = (k r ,k < j ) ) (e.g. 
Binney & Tremaine 2008, §6.2.2). Moreover, the compu¬ 
tation of the dynamical response of particles within our 
disc to a tightly-wound perturbing potential that oscil¬ 
lates at angular frequency v is covered by standard texts 
(see, for example, Binney & Tremaine 2008 §6.2.2(d) or 
Binney 2013 §4.2 for two different approaches). The re¬ 
sult is that a spiral potential ip^ (x)e 11 '* creates a spi¬ 
ral perturbation in the surface density of test particles, 
which through Poisson’s equation creates a response po¬ 
tential %p'(x)e wt that differs from the original stimulating 
potential only in magnitude. In fact 

=A k V> } (!?,</>), (30) 


where 


Ak = 


27rGS|fc r | 
k 2 (1 — s 2 ) 


H*,x)- 


(31) 


Such a packet is a non-trivial superposition of logarith¬ 
mic spirals, so M(v) cannot be diagonal in the basis pro¬ 
vided by logarithmic spirals. Physically, the dynamics of 
the disc is inherently local on account of the existence 
of resonant radii, so basis functions such as logarithmic 
spirals that extend from the disc’s centre to infinity can¬ 
not make M(v) diagonal. Given that we want M(y) to 
become diagonal, we must work with basis functions ip ^ 
that are local. 

Fouvry et al. (2014a) show how to construct a biothorg- 
onal basis of localised spirals. They divide the range 
(Rmin, Umax) of relevant radii into intervals of width 
a centred on R 0 , and then for any given wavevector 
k = (fe r ,fe^) create a basis function for each interval. 
Specifically, their basis potentials are 


V> (k ’ flo) (i?, </>) 


G 


\K\R 0 


gi {k r R+k^tj>) 
(7T<7 2 )V4 


exp 


(i?-i?o) 2 ] 

2<t 2 


(35) 

with k r Ro ^$> 1. The corresponding surface densities are 
E (k,flo) = (36) 


Here E is the disc’s surface density, 

K 


They show that two basis functions i /A k ’ fl o) and f/A k ,fl o) 
will be biorthogonal only when ARq = Rq — Rq and 
(32) A&y = kp — k 2 satisfy 


is the ratio of the frequency at which a star experiences 
the perturbation to the epicycle frequency, and T is the 
reduction factor (Kalnajs 1965; Lin & Shu 1966) 


Hs,x) = 2(1 



Xm r (x) 

1 — s 2 /m 2 ’ 


(33) 


where X mr is a modified Bessel function and the dimen¬ 
sionless quantity 


X = 


of k 2 

K 2 


(34) 


is a measure of how warm the disc is. In cases of interest 
the reduction factor is a number slightly smaller than 
unity and of little interest. 

The proportionality (30) suggests that in a basis 
formed of tightly-wound spiral waves the Fourier trans¬ 
formed response operator M(y) is diagonal with Ak the 
diagonal element associated with the given spiral wave. 
Hence the natural procedure might seem to be the adop¬ 
tion of the complete set of logarithmic spirals (e.g. Bin¬ 
ney & Tremaine 2008, §2.6.3) as the basis ip( p \ Un¬ 
fortunately, M(y) is not, in fact, diagonal in the basis 
formed by logarithmic spirals for the following reason. 
The demonstration that a spiral perturbation generates 
a spiral response scaled by Ak is a local result: the disc 
is analysed in just an annulus, and in the spirit of WKB 
analysis, the wave considered is a packet of finite length. 
Since the frequencies k and Cl# that appear in equation 
(32) are functions of radius, Ak is also a function of ra¬ 
dius whereas a diagonal element of M(y) should be a 
constant. Hence only a short packet of spiral waves pro¬ 
vides a good approximation to an eigenfunction of M(y). 


( AR 0 » a , or A R 0 = 0 , 

< 1 (37) 

| A k r — , or A k r = 0 . 

That is, the centres of neighbouring bands have to be sep¬ 
arated by more than the width of a band, and within any 
band, adjacent wavenumbers have to differ by enough to 
give a significant phase difference across the band. 

Now that the basis potentials have been chosen, one 
can use the mapping (x, v) i— > ( 0 , 3 ) provided by the 
epicycle approximation (e.g. Binney 2013, eq. 82) to com¬ 
pute their Fourier transforms with respect to the angle 
variables: 


V’, 


(k.Jto) 


(J) = drn+e R 


G 


A k r R s 


\kr\R 0 (7rcr 2 ) 1 / 4 


X Jm r ( V ^k, 


exp 


(Rg Ho)' 


2 a 2 


(38) 

where Rg(J<f>) is the radius of the circular orbit with an¬ 
gular momentum J$ and J mr is a Bessel function of the 
first kind. On account of the tight-winding condition 
k r R g 1, the phase shift 0° R is given by 


9r - -tt/2 • 


(39) 


In this basis, the response matrix M is diagonal, having 
diagonal elements 


-^(ki,fli)(k 2 ,K§) “ Kl $kt S nl ^ k ’ 


(40) 


where Ak is defined by equation (31). This expression for 
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M allows us to rewrite equation (29) in the form 


so by equation (28) 


i(J 


■o =Etfwi’‘[ ] ) ( i _ Ak ,) - (4i) 


In the chosen basis the expansion coefficients of the 
stimulating field are 


a p (v) = - J d 2 x[£ (p) (x)] *ip e (x, v) 


1 


Y GR P (tR7 2 ) 1/4 

x exp 


(R-R p 0 ) 


d RR 

21 


(42) 


2a 2 


e lRk -ip e kP (R,u), 


where 


r k p(R, = d <t>e~ ik >r{R, <j>, v) (43) 


is the Fourier transform in azimuthal angle and time of 
the stimulating potential. We further define the local 
radial Fourier transform of ij) e within the segment centred 
on by (Gabor 1946) 


A P q{v) 


_ l^vl^O (2 tt)- ik P( R P_ R 9) 

G (na 2 ) 1 / 2 

x 5(k p - k q )C(k p , v, (R p + R q 0 )/ 2). 


(49) 


To obtain the diffusion coefficients from equations (7) 
we substitute into equation (41) for c m our expressions 

(38) for ^ (J) and the above expression for A. We have 


Cm(J, V) = ^—y- E 


„k« 


pq 


k p) e -i k p AR p o-R q o) 


x exp 


(i? g - R p ) 2 + (R g - R q ) 2 


x 5(k p - A 9 ) 


2 a 2 

C(k p ,v, (R p + R q )/2) 


(l-A kP ) 2 


(50) 


The sums over p and q expand into sums over k p , k 9 , 

k p k q 

R p and Rq. On account of the factors 5^ and 5m,\ the 
sums over the (f> components of the k 1 are trivial. The 
sums over the k‘ r and R l 0 we approximate with integrals 
by the substitutions 


V’k* ( R o, v ) J dR exp 

x e- iiR - R ° )k r&p (R, v). 

**<f> 

This definition is motivated by the consequence that thus 
defined the local radial Fourier transform of a uniform 
potential ip e = 1 is independent of Rq. If in equation (42) 
we approximate the leading factor R in the integrand by 
R 0 , we then have 


(R-R p 0 Y 


2 cr 2 


(44) 


= 


Wl 


27T 


„-i R p n kl„ 


G (7ra 2 ) 1 / 4 


o k H’UK")- ( 45 ) 


We require the ensemble average a p (z/)a*(i/) (eqs. 28 
and 29), which is related to the ensemble average 
^> e (x, f)-i/> e (x', t'). We assume that stimulating fluctua¬ 
tions are quasi-stationary in the sense that 

r k 2 R A)r k *2R'A') = c k ,(t -t^R-m, (r+ my 2 ), 

(46) 

with the dependence on R + R' being weak. With this 
assumption that the process ip^ is stationary in time 
and “locally stationary” in space, it follows that 

$u.p( r , v)$ kq {R', v') = 2n5(v - v')5(k p - A 9 ) ^ 

x C'(k p ,i/,(i? + i?')/2), 


where C (k, u, R) is the spatio-temporal power spectrum 
of the stimulating noise in the neighbourhood of R. Now 
we can write 


a P (v)a* q (y') 


_ \K\K (2?r) 3 j K(R p n-R q „) 5(l . _ u >\ 

~ G (na 2 ) 1 / 2 1 ’ 

x5(k p -k q )C( k*>, (R p + R q )/ 2), 

(48) 


E/(M 

kqp 

Es( fi o) 

Ro 


1 

A Ay 

1 

A 


dk r /(Ay), 


di?o g(Ro), 


(51) 


where Ak r is the difference between successive values of 
Ay in the sum, and similarly for ARq. Then the Dirac 
delta function in equation (50) allows us to integrate over 
A 9 . We assume that cr is small enough that we can ap¬ 
proximate each Gaussian exponential by \/2na5(R g —R} ) ) 
so we can trivially integrate over the R l 0 . Finally, the 
presence in equation (48) of a rapidly oscillating com¬ 
plex exponential e lkrR ° imposes that the intervals AAy 
and Ai? 0 must satisfy the critical-sampling condition 
Ak l r ARl = 27t (Daubechies 1990). With this condition, 
equation (50) simplifies to 

f /Jm r k2 

c m (J, (c) = 2 J dk r - (1 _ Ak)2 7 C(Kv,R s ), (52) 

where by hypothesis the dependence of C on R g is weak, 
and we have k<f, = my This is our principal result. It 
enables us to compute the diffusion tensor at any point in 
the action space of a thin, self-gravitating disc given the 
power spectrum of the noise that excites spiral structure 
in the disc. 


4. APPLICATION TO A MESTEL DISC 

We apply our results to the same Mestel disc (Mes- 
tel 1963) that S12 discussed. The basic properties of 
this disc are given in §2.6.1(a) and §4.5.1 of Binney & 
Tremaine (2008). Its circular speed is a constant Vo and 
its potential is 


i>{R) = Vo ln 


R 

Rin 


(53) 
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where the value of R[ n is arbitrary. The corresponding 
surface density is 


£ (R) 


v 2 

v o 

2nGR ' 


(54) 


Toomre (1977) gives DF f(E, J^) that self-consistently 
generates this disc. When we use the epicycle approxi¬ 
mation to replace the energy E by the radial action J r , 
the DF becomes 


fo(Jr, J*) = G9(J 0 ) 
where 
C = 


Jd> 


RinV 0 


exp - 


k(J<p ) 


(55) 


V 0 q+2 /(2rrGR in ) 


and 


2®/ 2 v / ^(g/2- l/2)\cr q+2 




exp(-V r 0 2 /2cr 2 ) , (56) 




(57) 


and O ( J^) is an Heaviside function removing retrograde 
stars. 

Since the central singularity and infinite extent of the 
Mestel disc are problematic, it is customary to modify 
the DF (55) by multiplying it by factors T( J^) that taper 
the stellar distribution at very small and very large radii. 
These factors are 


Tin{J<l>) 


T ou t(J<p ) 


(RinVoY + j; ’ 



(58) 


where v and fi control the sharpness of the two tapers. 
T; n models the presence of a bulge by diminishing the 
DF inward of R ln . Here, T out models the outer edge of 
the disc, beyond which the gravitational field is entirely 
generated by dark matter. Even after tapering the stel¬ 
lar distribution, ip(R) continues to be given by equation 
(53) because the bulge and the dark halo are presumed 
to provide the gravitational force that was originally pro¬ 
vided by the un-tapered disc. 

In our numerical work we use the same taper constants 
as S12. We adopt a system of units such that: Vo = 
G = R[ n = 1. The other numerical factors are given by 
q = 11.4, v = 4, /i = 5, Tout = H-5. 

Within the epicyclic approximation, the azimuthal and 
radial frequencies are 


fi W = -Jr ’ 

Kg 

k(J</>) = v^fi(J^), 


(59) 


and are thus independent of J r . The ratio k/Q = \[2 
is a constant. This ratio determines the location of the 
resonances, so it is important for the disc’s dynamical 
behaviour. By taking it to be a constant we risk intro¬ 
ducing unphysical artifacts in the dynamics. 

The distribution function (55) multiplied by the taper 
factors (58) takes the form of a locally isothermal -DF or 





Figure 1. Surface density Et of the tapered Mestel disc. The 
unit system has been chosen so that Vo = G = = 1. 



0 2 4 6 8 


J(!> 

Figure 2. Contours of the initial distribution function in 
action-space (J^>,J r ), within the epicyclic approximation. 
The contours are spaced linearly between 95% and 5% of the 
distribution function maximum. 


Scliwarzschild-DF with the correct normalization, which 
can be rewritten as 


Wr, J<f>) 


fi (J»)St(^) 

n k(J^) a 2 1 



(60) 


where the intrinsic frequencies are given by equation (59) 
and the taped surface density in analogy with equa¬ 
tion (54) is given by 

St(j0) = 2^au 9(J0) Tin(j0) ToM) ' (61) 


The shape of the damped surface density is shown in 
Figure 1. Figure 2 shows the level contours of the distri¬ 
bution function fo- 

In the scale-invariant Mestel disc the local Toomre 
(1964) parameter 


(Jr fi:(T<^) 

3.36 GE^) 


(62) 


is independent of radius, and in the tapered disc Q is 
correspondingly flat between the tapers. For realistically 
small values of oy/Vo, the plateau in Q can lie below 
unity, making the disc unstable. To keep Q everywhere 
well above unity it is conventional to suppose that only 
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Figure 3. Variation of the Toomre parameter Q with the 
angular momentum J$. It is scale invariant except in the 
inner/outer regions because of the presence of the tapering 
functions Tin and Tout- The unit system has been chosen so 
that Vo = G = Rin = 1- 


a fraction £ < 1 of the disc is self-gravitating with the 
rest of the gravitational field provided by an unresponsive 
halo. In the S12 simulation, the fraction of active surface 
density was £ = 0.5. The dependence of Q with radius 
with this value of £ is shown in Figure 3 - Q ~ 1.5 
between the tapers and increases strongly in the tapered 
regions. 


4.1. Impact of shot noise 

To proceed further we need to assume some form for 
the power spectrum C{ k, v , R g ) that appears in equation 
(52). An inevitable source of noise is shot noise caused 
by the the finite number of stars in the disc, and mas¬ 
sive, compact gas clouds are a source of spectrally similar 
noise, so let us investigate the impact that shot noise has. 
In this case the power spectrum is independent of v and 
k r , and varies with radius like ^/E(i?). Then C oc E, 
so to within a normalisation that depends on particle 
number, we have 

__ /' 3m r (\/2Jr/K kr\ 

Cm(J, f) = Et(J^) J d k r - ^ -. (63) 

The eigenvalues Ak have to be evaluated at v = mfl, and 
then s = m r by equation (32). In order to handle the 
singularity of the equation (31) when s = ±1, one adds 
a small imaginary part to the frequency of evaluation, 
so that s = in r + i rj. As long as 77 in modulus is small 
compared to imaginary part of the least damped mode 
of the disc, adding this complex part makes a negligible 
contribution on the expression of Re(A). 

In equation (63) the integral over k r should formally 
be over the full range 0 to oo. However, small values of 
k r are unphysical and violate our assumption of tightly- 
wound spirals. Values of k r that are larger than ~ 2tt 
divided by the thickness of a galactic disc are also un¬ 
physical, and in the case m r = 0 of the CR the integral 
diverges at J r = 0 since then the Bessel function remains 
non-zero to arbitrarily high k r and (1 — Ak) -2 is always 
greater than unity. Hence we must determine appropri¬ 
ate upper and lower limits to the integration on k r . 

At any point in action space the biggest contribution to 
the diffusion tensor will come from waves that yield the 


A 



Figure 4. Variation of eigenvalues A of the response matrix 
with the WKB-frequency k r for two values of J^. The curve 
that peaks at small k r is for the larger value of J^>. 


largest value of Ak- Hence we now examine the structure 
of the function k r >->• Ak for given m and J. Figure 4 
shows that Ak has a well-defined peak at a value k max of 
k r that decreases as J$ increases. In fact /c max ex 1 / 
because the radius of a near-circular orbit is R oc and 
in a scale-free model we expect fc max oc 1 /R. For the 
same reason we expect the width A*, of the peak in Ak 
to be proportional to 1/ In light of these observations 
we adopt as lower and upper bounds on the k r integral 
the wavenumbers k ln f and fc sup defined by 

Afci„f = A fcsup = ^ A fcmax . (64) 

These limiting values of k r are marked on Fig. 4. 

At each point in action space there are contributions 
to the diffusion coefficients from several values of m. 
The contribution from m r = — 1 is driven by waves 
that have their inner Lindblad resonance at that point, 
while that from m r = +1 is driven by waves that have 
their outer Lindblad resonance there, and the contribu¬ 
tion from m r = 0 is driven by waves locally in corotation. 
Since Ak depends on s 2 = to 2 , the value of Ak is the same 
for m r = ±1. At a given point in action space different 
values of the frequency v are associated with m r = ±1, 
but since we are considering shot noise, the fluctuations 
have the same power at all frequencies. Hence the values 
m r = ±1 contribute equally to the diffusion coefficients. 
The lower curve in Fig. 5 shows the extent to which this 
contribution is amplified by the disc’s self-gravity and we 
see that the amplification is modest. The upper curve in 
Fig. 5 shows the amplification by self-gravity in the case 
m r = 0 of corotation: it is much larger. 

In a cool disc, \df /dJ r \ \df /dJ^\ so by equation (5) 
a reasonable approximation to the diffusive flux, when 
m r ^ 0, is 

Fj = - \ TO,;C m TO r ^y-. (65) 


It follows that waves at inner Lindblad resonance, which 
couple through (m r ,m ( f,) = (—1,2) drive a diffusive flux 


F 


_ ( _ 1> 2 )5 C (~1,2) 


dfo 

dJ r 


( 66 ) 


while waves at outer Lindblad resonance drive the flux 


F — (1) < F)\ c (i,2) 


dfo 

dJ r 


(67) 
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1/(1 -A max ) 



Figure 5. Dependence of the amplification factor 1/(1 — 
Amax) with the position Throughout the disc, the am¬ 
plification of waves at corotation is larger than that of waves 
at inner Lindblad resonance. 

These formulae show that waves with a Lindblad res¬ 
onance at J drive a flux towards increasing J r : these 
waves are heating the disc by increasing the eccentrici¬ 
ties of orbits. Wave that have their ILR at J drive stars 
to lower angular momentum, while those with their OLR 
at J drive stars to higher angular momentum. We shall 
see below that at low waves with a local ILR domi¬ 
nate, so on average angular momenta decrease, while at 
large J the waves with a local OLR dominate and angu¬ 
lar momenta increase on average. Thus spiral structure 
conveys angular momentum outwards, thus liberating en¬ 
ergy that is used to heat the disc. 

Waves at corotation couple through m = (0, 2) so (a) 
they can only drive diffusion parallel to the axis, and 
(b) that diffusion is proportional to the gradient’s small 
component df/dJ In practice this diffusion is impor¬ 
tant only in so far as the disc has a metallicity gradient, 
when the DF of stars of any particular metallicity can 
have a significant derivative with respect to and the 
large value of c m at corotation can have a big impact on 
the metallicity distribution - radial migration at coro¬ 
tation is important for chemical evolution (Sellwood & 
Binney 2002; Schonrich & Binney 2009). 

4.2. Reproducing the S12 simulation 

In this section we examine the extent to which our 
analytic results can explain the simulations of tapered 
Mestel discs described by S12. We do not expect perfect 
agreement between our results and the numerical experi¬ 
ments because we have employed several approximations. 
In particular we have assumed that the driving fluctua¬ 
tions are white noise that has an amplitude squared that 
is proportional to the disc’s surface density. We have 
assumed, moreover, that these fluctuations drive tightly- 
wound waves. S12 restricted disturbing forces to = 2, 
so we impose this same restriction on m. 

The dark contours in Fig. 6 show the magnitude of the 
diffusive flux that is generated by the two Lindblad res¬ 
onances and corotation when one adopts the same nu¬ 
merical pre-factors as S12. The gray arrow shows the 
direction of the diffusive flux at the location of the peak 
flux; the direction is similar at neighbouring points. The 
thin contours in Fig. 6 show the value of the distribution 
function in S12 after diffusion has taken place. Originally 
the DF peaked along the axis at ~ 1, becoming 


J R 


0.4 



Figure 6. Map of the norm of the total flux summed over 
the three resonances (ILR, CR, OLR) (bold lines). The con¬ 
tours are spaced linearly between 95% and 5% of the function 
maximum. The gray vector gives the direction of the vector 
flux associated to the norm maximum (arbitrary length). The 
background contours correspond to the diffused distribution 
from S12 (thin lines), which exhibits a narrow ridge of diffu¬ 
sion. 

small to the left of that point on account of the inner 
taper. Above the location of the peak DF, there is a 
prominent horn of enhanced stellar density that extends 
roughly in the direction of the diffusive flux. Thus our 
analytic results are in good qualitative agreement with 
the numerical experiments of S12. 

In Fig. 6 waves at ILR drive diffusion parallel to vectors 
that are inclined by 153° to the axis, whereas the net 
flux makes an angle of 111° with the axis. That these 
angles are similar attests to the dominant role of waves 
that have ILR near where the DF peaks. This dominance 
arises because dfo/dJ r is always negative and in this 
region dfo/dJ$ > 0, so |m • dfo/dJ\ is larger for the 
waves at ILR, which have m r = — 1, than for the waves 
at OLR, while c m is the same for both values of m r . 

Given that we have assumed that the driving fluctu¬ 
ations are white noise, it is perhaps surprising that the 
magnitude of the diffusive flux is as sharply peaked in 
action space as the dark contours in Fig. 6 show it to 
be. This localisation surely reflects the sharp tapering of 
the DF at small radii. Just outside this taper the disc 
becomes significantly self-gravitating and most effective 
at amplifying the stimulating white noise. The amplified 
waves tend to have their ILRs further in, where the taper 
is cutting into the disc. 

We have seen that self-gravity amplifies stimuli most 
strongly at corotation (Fig. 5). In fact, the dominance of 
corotation increases without limit as Q decreases because 
the value of A max associated with corotation approaches 
unity, and the associated amplification diverges, while 
the value of A max associated with the LRs remains sig¬ 
nificantly below unity. Fig. 7 illustrates this phenomenon 
by comparing the velocity of diffusion in action space, 


F 



(68) 
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for two values of For the value £ = 0.5 assumed by S12 
(left panel) the velocity vectors are significantly non-zero 
only within the narrow strip at ~ 1 associated with 
the ILR and point up and to the left. For £ = 0.73 the 
velocity is also quite large at > 2 near the axis and 
is there directed up and to the right. 

5. CONCLUSIONS 

In a star cluster or a galaxy stars move on orbits in 
the system’s mean gravitational field, but the orbit each 
star is on evolves slowly in response to fluctuations in 
the mean field (Weinberg 2001a). In a hot stellar sys¬ 
tem such as a star cluster, the dominant fluctuations are 
pure shot noise arising from the finite number of stars in 
the system. In a disc galaxy the situation is much more 
complex and interesting because the system’s response is 
more frequency-sensitive, and as Toomre’s Q —> 1 stim¬ 
uli are strongly amplified and distorted by the coherent 
responses they induce in the disc. 

The orbit-averaged Fokker-Planck equation, which de¬ 
scribes the evolution of the distribution function as stars 
diffuse through action space, provides the mathematical 
device of choice to compute the long-term evolution of a 
stellar system. The equation is readily solved once the 
diffusion tensor has been computed. 

We have laid out the general formalism for computing 
the diffusion tensor in the presence of significant coherent 
response to stimuli, and have implemented this formal¬ 
ism in the case of a razor-thin, cool disc. By introduc¬ 
ing a set of basis functions for the disc’s potential that 
comprise sets of localised, tightly wound spirals, we have 
derived equation (52), which reduces computation of the 
diffusion tensor to execution of a one-dimensional inte¬ 
gral over the auto-correlation function of the stimulating 
noise. 

We used this equation to compute the diffusion tensor 
for a Mestel disc that has been trimmed at both large and 
small radii by tapers and is exposed to shot noise. We 
find that diffusive flux shows quite sharp peaks in action 
space. When Toomre’s parameter Q is significantly big¬ 
ger than unity, the diffusive flux is quite localised near 
the inner taper, and is primarily driven by waves that 
have their inner Lindblad resonances there. This result 
is interesting in relation to the N-body simulations of 
S12 because in these simulations, orbital diffusion led to 
the formation of a horn of enhanced density in phase 
space that lies close to the region in which our analytic 
work predicts that the diffusive flux peaks. Moreover, 
the direction along the horn is similar to the predicted 
direction of the diffusive flux. Thus it appears that our 
analytic work recovers quite well the main feature of the 
simulations. 

As Q approaches unity, the corotation resonances of 
waves become more important because self-gravity am¬ 
plifies perturbations at corotation very strongly. Waves 
that are in corotation at some point in action space drive 
diffusion parallel to the angular-momentum axis of ac¬ 
tion space, i.e. they drive radial migration. By contrast, 
waves that are at a Lindblad resonance heat the disc by 
driving stars to higher eccentricity, while either decreas¬ 
ing (at ILR) or increasing (at OLR) angular momentum. 
Consequently, our calculations predict that the amount 
of radial migration for a given level of heating increases 
as Q decreases towards unity. 


In our application of action-space diffusion to cool 
discs, the disc’s response to stimuli has been restricted to 
tightly-wound spirals. This restriction unfortunately ex¬ 
cludes from consideration the key physics of “swing am¬ 
plification” at corotation (Toornre 1981). As Goldreich 
& Lynden-Bell (1965) originally showed in the context of 
a gas disc, leading waves are amplified at corotation as 
they morph into trailing waves. These waves will then 
propagate to the Lindblad resonances and there heat the 
disc (Toomre 1981). Shot noise will stimulate leading 
waves in the same amount as trailing waves, and the 
leading waves will move to corotation rather than to the 
Lindblad resonances, and there morph into trailing waves 
of larger amplitude. Our computation has not included 
the effect of swing-amplification at corotation, and will 
consequently under-estimate the extent of heating. The 
under-estimation will be largest for small values of Q , 
because swing amplification diverges as Q —> 1. 

Another shortcoming of the present analysis is the 
crude noise model which has ad-hoc dependence on the 
temporal frequency v and the radial wavenumber k r , and 
is only a function of the position in the disc through 
The model aims to reproduce the Poisson shot noise 
caused by the finite number of particles in the disc. In a 
companion paper, we will investigate the WKB limit of 
the Balescu-Lenard equation. This approach will allow 
us to avoid such approximation since it naturally cap¬ 
tures the intrinsic noise, due to finite— N effects, and its 
impact on the quasi-stationnary distribution function, as 
long as the evolution of the system is made through tran¬ 
sient tightly wound spirals. 

In a forthcoming paper we will use the present formal¬ 
ism of forced diffusion to study the evolution of a stellar 
disc as a result of cosmic noise: the noise generated by 
satellites that orbit in a disc’s host dark halo. 

Another possible application of the formalism devel¬ 
oped here is to study the secular evolution of dark-matter 
cusps at the centres of galaxies in response to stochastic 
excitation by the inner baryonic disc and bulge. (Fouvry 
et al. 2014b). 

This paper is dedicated to the memory of Jean Heyvaerts 
who was its inspiration. 
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